#include"stationary_solver.h"

StationarySolver* StationarySolver_ = NULL;

StationarySolver::StationarySolver(GroundData* GroundData_, PipeData* PipeData_) {
  this->GroundData_ = GroundData_;
  this->PipeData_ = PipeData_;
  T_new = new float*[GroundData_->cell_q_x];
  for(int i=0; i<GroundData_->cell_q_x; i++)
    T_new[i] = new float[GroundData_->cell_q_y]; 
  Q = new float[GroundData_->cell_q_y];
  c = 0.003f;
  k = 1.0f;
  ro = 0.01f;
  return;
}

int StationarySolver::SetGroundInitialConditions(float T_0, float T_0_surface, float T_0_earth) {
  this->T_0_earth = T_0_earth;
  this->T_0_surface = T_0_surface;
  for(int i=0; i<GroundData_->cell_q_x; i++) {
    for(int j=0; j<GroundData_->cell_q_y; j++) {
      GroundData_->T[i][j] = T_0;
      //GroundData_->T[i][j] = T_0+5-i;
      if(j == 0)
        GroundData_->T[i][j] = T_0_surface;
    }
  }
  return 1;
}

int StationarySolver::SetPipeInitialConditions(float T_in, float T_out, float v) {
  for(int i=0; i<PipeData_->cell_q_y; i++) {
    PipeData_->T_outside[i] = T_in + i*(T_out-T_in)/(2.0f*PipeData_->cell_q_y);
  }
  for(int i=1; i<=PipeData_->cell_q_y; i++) {
    PipeData_->T_inside[PipeData_->cell_q_y - i] = T_in + (i+PipeData_->cell_q_y)*(T_out-T_in)/(2.0f*PipeData_->cell_q_y);
  }
  this->v = v;
  return 1;
}

int StationarySolver::MakeStep(float delta_t) {
  CountStream();
  float r_now = 0.0f;
  // main part
  for(int i=1; i<(GroundData_->cell_q_x-1); i++) {
    r_now += GroundData_->delta_x[i];
    for(int j=1; j<(GroundData_->cell_q_y-1); j++) {
      T_new[i][j] = GroundData_->T[i][j] + 
        delta_t*k*(/*d2t*/(GroundData_->T[i-1][j] - 2*GroundData_->T[i][j] + GroundData_->T[i+1][j])/(GroundData_->delta_x[i]*GroundData_->delta_x[i]) +
        1.0f/r_now * (GroundData_->T[i+1][j] - GroundData_->T[i][j])/GroundData_->delta_x[i] + (GroundData_->T[i][j-1] - 2*GroundData_->T[i][j] +
        GroundData_->T[i][j+1])/(GroundData_->delta_y*GroundData_->delta_y));
    }
  }
  // stream x=0
  r_now = 0.0f;
  for(int i=0; i<1; i++) {
    r_now += GroundData_->delta_x[i];
    for(int j=1; j<(GroundData_->cell_q_y-1); j++) {
      T_new[i][j] = GroundData_->T[i][j] + 
        delta_t*k*(/*d2t*/(GroundData_->T[i][j] - 2*GroundData_->T[i][j] + GroundData_->T[i+1][j])/(GroundData_->delta_x[i]*GroundData_->delta_x[i]) +
        1.0f/r_now * (GroundData_->T[i+1][j] - GroundData_->T[i][j])/GroundData_->delta_x[i] + (GroundData_->T[i][j-1] - 2*GroundData_->T[i][j] +
        GroundData_->T[i][j+1])/(GroundData_->delta_y*GroundData_->delta_y));
      T_new[i][j] += delta_t*Q[j]/(c*ro);
    }
  }
  // x=cell_q_y-1
  for(int j=1; j<(GroundData_->cell_q_y-1); j++) {
    T_new[GroundData_->cell_q_x-1][j] = GroundData_->T[GroundData_->cell_q_x-1][j];
  }
  // ground surface (up)
  r_now = 0.0f;
  for(int i=0; i<GroundData_->cell_q_x; i++) {
    r_now += GroundData_->delta_x[i];
    for(int j=0; j<1; j++) {
      T_new[i][j] = T_0_surface;
    }
  }
  // ground on depth
  r_now = 0.0f;
  for(int i=0; i<GroundData_->cell_q_x; i++) {
    r_now += GroundData_->delta_x[i];
    for(int j=GroundData_->cell_q_y-1; j<GroundData_->cell_q_y; j++) {
      T_new[i][j] = GroundData_->T[i][j];
    }
  }
  // reload T[][]
  for(int i=0; i<GroundData_->cell_q_x; i++) {
    for(int j=0; j<GroundData_->cell_q_y; j++) {
      GroundData_->T[i][j] = T_new[i][j];
    }
  }
  return 1;
}

int StationarySolver::CountStream(void) {
  for(int j=0; j<PipeData_->cell_q_y; j++) {
    Q[j] = c*v/PipeData_->delta_y*(PipeData_->T_outside[j] - PipeData_->T_outside[j+1] + 
      PipeData_->T_inside[j+1] - PipeData_->T_inside[j]);
  }
  return 1;
}